Isogeometric Methods for Computational Electromagnetics: 
B-spline and T-spline discretizations 



A. Buffa a , G. Sangalli b ' a , R. Vazquez a 

a Istituto di Matematica Applicata e Tecnologie Informatiche 'E. Magenes' del CNR 

via Ferrata 1, 27100, Pavia, Italy 
b Dipartimento di Matematica, Universita di Pavia, via Ferrata 1, 27100, Pavia, Italy 



Abstract 

In this paper we introduce methods for electromagnetic wave propagation, based on splines and on 
T-splines. We define spline spaces which form a De Rham complex and, following the isogeometric 
paradigm, we map them on domains which are (piecewise) spline or NURBS geometries. We 
analyse their geometric structure, as related to the connectivity of the underlying mesh, and we 
give a physical interpretation of the fields degrees-of-freedom through the concept of control fields. 
The theory is then extended to the case of meshes with T-junctions, leveraging on the recent theory 
of T-splines. The use of T-splines enhance our spline methods with local refinement capability and 
numerical tests show the efficiency and the accuracy of the techniques we propose. 

Keywords: Maxwell equations, De Rham diagram, Exact Sequences, Isogeometric Methods, 
Splines, T-splines. 



1. Introduction 

Electromagnetic field computations and, more generally, the numerical discretization of equa- 
tions enjoying a relevant geometric structure, is one of the most interesting challenge of numerical 
analysis for PDEs and several results have been obtained in the last decade. Indeed, only for 
Galerkin methods, three Acta Numerica overview papers have been published: by Hipmair [1] , by 
Arnold, Falk and Winther [2J, and by Boffi [3], addressing different aspects of the problem. 

On the one hand, discrete schemes have to preserve the geometric structure of the underlying 
PDEs in order to avoid spurious behaviors, instability or non-physical solutions (see e.g., the 
pioneering paper jl]). For electromagnetics, as it is clear from the references above, numerical 
schemes have to be related with a discrete De Rham complex. On the other hand, especially in 
view of high frequency computations, numerical schemes need to be efficient and accurate. This 
requires many features, and among others it requires adaptivity, or at least local mesh refinement 
capability, in order to capture the strong singularities of the electromagnetic field, possibly driven 
by a-posteriori error estimator as, e.g., in [5]. 

In this paper we present and analyse discretization techniques for electromagnetic fields based 
on splines and generalizations of splines, as NURBS ([6]) or T-splines ([7] or below). Our work 
originates from IsoGeometric Analysis (IGA), [8J. Isogeometric analysis has been introduced in 
2005 by Hughes and co-authors in the seminal paper jH] to solve structural mechanic problems 
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directly on the geometry output by a CAD system, and has set the paradigm to use splines, 
NURBS or their generalization as generating functions for the construction of Galerkin spaces. 
This idea has been proved to be extremely effective and IGA is spreading very fast across different 
scientific communities: structural mechanics (see e.g., [8J, [10], [II] . [12], [13], [H], [15]), geometric 
modeling (see e.g., [IB] . [17] , [IK] and also [TU]) and numerical analysis (see e.g., [20], [21] , [22], 

[23]. m, m, m.)- 

In this paper we present the recent advances in the use of the isogeometric paradigm and 
spline-based methods for electromagnetic wave computations. This research has started with the 
two papers [2Z] and [22] and can likely be considered as still in infancy (see also [2H] for the 
applications of this results). This paper aims at showing the potential impact of these techniques 
in the electromagnetic community by addressing several aspects: from the geometric structure of 
the proposed methods, to local refinement strategies. 

We introduce the spline complex studied in [22] (see (38) and (39)) and we present its 
properties: we construct canonical bases so that the matrices representing differential operators 
are the incidence matrices of the underlying meshes, and this enlightens the relation between the 
spline complex and the geometry of the underlying meshes. We show that for different choices of 
the degree of splines, the spline complex is isomorphic to the co-chain complex or to the chain 
complex of the underlying mesh. Besides this interesting fact, we also introduce the concept of 
control fields in analogy to control points which are ubiquitous in spline theory (see e.g., [21] or 
[30] ) which provide the correct physical interpretation of degrees of freedom. Finally, we extend the 
results of [22] to multi-patch geometries, i.e., geometries which are piecewise spline or NURBS 
mapping of the unit cube. We refer the reader to [26] for a detailed description of this class of 
geometries. 

The second major contribution in this paper is a step towards adaptivity for spline-based 
methods. Leveraging on the recent work on T-splines, we design a two dimensional T-spline 
complex where meshes with T-junctions can be used to allow for adaptivity. T-splines are the 
most attractive way to break the tensor product structure of splines while keeping their structure 
and their accuracy. T-splines have been introduced in [7] and [31] and their use as a fundamental 
tool to enhance isogeometric analysis with adaptivity has been proposed in [32] . A series of papers 
has followed [33], [31], [35], [36], together with the relevant class of Analysis Suitable T-splines 
[3"7] . [55] which we use in our construction. The two dimensional T-spline complex is used to 
treat three dimensional problems with symmetry. We should also mention that the definition and 
use of T-splines in three dimension are not yet well understood, but object of an intensive study. 
Their use will allow, on a longer time perspective, to design full adaptive algorithms, on very 
general geometries parametrized on totally unstructured meshes. We refer the reader to [39] for a 
monograph on the modern use of T-splines in geometric modeling. 

Finally, we should remark that the spline spaces we study in this paper have a wide domain 
of applications and can be applied successfully to the discretization of other problems than elec- 
tromagnetics. In fact, they can be used to solve the Darcy flows equations or more generally the 
Hodge laplacian operator as detailed in [2] and [ID]. Moreover, thanks to the regularity of spline 
spaces, their use in fluids is very promising. In the paper [UJ they are used for the first time to 
solve the Stokes equations, in [42] the Stokes eigenvalue problem is addressed, and in the sequence 
of three papers [13], [B], jl5] they are applied to solve Stokes and Brinkman equations, steady 
and unsteady Navier-Stokes equations, providing impressive results. 

The outline of the paper is the following. In Section[2]we set up the notation for the problems we 
address, in Section[3]we present known results about splines and NURBS in a self-contained way; in 
Section [4] we present the spline complex and all the related results while in Section [5] we introduce 
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the T-spline complex and analyse its properties. Finally, in Section [6] we present numerical results: 
the first ones are two and three dimensional, academic tests aiming at demonstrating the validity 
of the proposed approach. As a last example, we compute the propagation in a waveguide with 
geometric inhomogeneity, on a three dimensional locally refined mesh. 



2. Notation 



In this section we present the notation that we need to describe the time-harmonic Maxwell 
problem. Let ft C IR 3 be a bounded Lipschitz domain. We denote by L 2 (ft) the space of complex 
square integrable functions on ft, endowed with standard L 2 norm || • ||l 2 (c), and by L 2 (ft) their 
vectorial counterparts. The Hilbert space H 1 ^) contains functions of L 2 (ft) such that their 
first order derivatives also belong to L 2 (ft). We denote by Hq(Q) C if 1 (ft) the subspace of 
functions with homogeneous boundary condition. We will also make use of the space H(curl; ft) , 
constituted by all functions in L 2 (ft) such that their curl also belongs to L 2 (ft), and H(div; ft) , 
the space of functions in L 2 (ft) such that their divergence belongs to L 2 (ft). Moreover, we denote 
by H (curl;ft) (resp. H (div; ft) ) the subspace of H(curl; ft) (resp. H(div; ft) ) of fields with 
vanishing tangential (resp. normal) component. 

For the sake of simplicity, we assume that the domain ft, referred to as physical domain in the 
following, is bounded Lipschitz and simply connected, and that its boundary <9ft is connected. We 
also assume that it is defined through a continuously different iable parametrization with contin- 
uously differentiable inverse which we denote as F : ft — > ft, where ft will be referred to as the 
parametric domain. Further assumptions on the geometrical mapping F will be given later. 

Some notation will be borrowed from the context of differential forms: first of all, we define 
the spaces 

X° := H\tt) , X 1 := H(curl; ft) , X 2 : = H(div; ft) , X 3 := L 2 (ft) , 
X° := # x (ft) , X 1 := H(curl; ft) , X 2 := H(div; ft) , X 3 := L 2 (ft) ; 

Since the parametrization F and its inverse are smooth, we can define the pullbacks that relate 
these spaces as (see [H Sect. 2.2]): 



, 2 (v) 

t %) 



0oF, 

(DF) T (uoF), 

det( J DF)( J DF)- 1 (voF) 

det(DF)(poF), 



0GX°, 
u G X 1 , 
vGX 2 , 

v?gx 3 , 



(1) 



where DF is the Jacobian matrix of the mapping F. Then, due to the curl and divergence 
conserving properties of l 1 and t 2 , respectively (see |HU Sect. 3.9], for instance), the following 
commuting De Rham diagram is satisfied (see [H Sect. 2.2]): 



-)• X c 



grad 



-)• x c 



grad 



> X 1 



curl 



> X 1 



curl 



> X 2 



div 



> X 2 



div 



> X 3 



> X 3 



-> 



(2) 



-+ 0. 



We are also interested in spaces with boundary conditions, denoted with the subindex 0, 
X ° := i^ft) , Xl := H (curl; ft) , X 2 := H (div; ft) , X 3 := L 2 (ft) , 



X °: = 



^(ft), X] := H (curl;ft), X 2 := H (div;ft), X 3 := L 2 (ft) 
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for which the De Rham diagram reads 
- 



o > x ° xl 



curl 



curl 



div 



div 



(3) 



which also expresses the integral preserving property of l 3 . 

Remark 2.1. As it is well known, the exactness of the sequences (J2J) and ^ relies on the assump- 
tion that Q (and Q) has a trivial topology. All what we develop in this paper applies in principle 
also to the case of arbitrary topology but we do not present all the details here. 



3. Preliminaries on splines and NURBS 

We give here a brief overview on B-splines and, in the spirit of [32], we also introduce some 
concepts that will be needed in the definition of T-splines. For more details on B-splines we refer 
the reader to [9j |6] . 

3.1. Univariate B-splines 

3.1.1. Knot vector and B- spline functions, refinement, spline derivatives 

Given two positive integers p and n, we say that H := {£i, . . . , £ n+p+ i} is a p-open knot vector 

if 

SI = • • • = £p+l < Sp+2 < • • • < in < £n+l = • • • = £ n+p+ l, 

where repeated knots are allowed and denote by rrij the multiplicity of the knot £j. We assume 
rrij < p + 1 for all internal knots. 

From the knot vector H, B-spline functions of degree p are defined following the well-known 
Cox-DeBoor recursive formula: we start with piecewise constants (p = 0): 

and for p > 1 the B-spline functions are defined by the recursion 

N ilP (0 = t^t^i(C) + /' +P+1 ~ C Wi(0- (5) 

Si+p — Si Si+p+1 — Si+1 

This gives a set of n B-splines that form a basis of the space of splines, that is, piecewise polynomials 
of degree p with p — rrij continuous derivatives at the internal knots £j, for j — p + 2, . . . , n. We 
denote this univariate spline space by 

S P (E) = span{iVi iP , i = 1, . . . , n} (6) 

An example of some B-splines is given in Figure [Tj Notice that the B-spline function Ni )P is 
supported in the interval an d in fact its definition only depends on the knots within 

that interval. For this reason, we define the local knot vector S» )P = {CijCi+i; • • • jCi+p+i}; an d we 
will sometimes denote N ijP (() = N[Ei tP ](Q. 

In the context of splines, three kinds of refinement are possible, as explained in [§]: 
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Figure 1: Example of B-splines of degree 2 (left) and 3 (right). 

1. k-refinement which corresponds to successive application of the Cox-DeBoor formula Q-Q. 
Regularity is raised together with the degree: therefore, the spaces ^ are not nested under 
/c-refmement but, at each step (degree and regularity elevation), the dimension of the space 
increases by 1. The name fc-refmement has been introduced in [9]; 

2. h-refinement which corresponds to mesh refinement and is obtained by knot insertion. Let 
H := {£i, . . . , £fc, £, 6c+i, £ n +p+i} be the knot vector after inserting the knot £ in S. Then, the 
new B-spline functions {Ni iP (Q, . . . , N n+p+2 , P (()} are constructed as follows: 

N itP (0 = a t N hp (() + (1 - a,)iV,_ liP (C) (7) 

where oij = 1, if 1 < i < k — p, ati — — - — ^— , if A; — p + 1 < i < k and, = for 

+ 1 < z < n + p + 2. When £ is equal to £fc or or to both, the knot insertion 
corresponds to reduction of the inter-element regularity at £. 

3. p-refinement which corresponds to the degree raising with fixed interelement regularity, and 
generates a sequence of nested spaces. 

Assuming the maximum multiplicity of the internal knots is less than or equal to p, i.e., the 
B-spline functions are at least continuous, the derivative of the B-spline N itP is a spline as well. 
Indeed, it belongs to the spline space 5 p _i(S'), where 5' = {£2, • • • ,£ n + P } is a (p — l)-open knot 
vector. Obviously, the regularity of splines in S p _i(E') is one less than the regularity in S P (E). 

In the following we assume that £1 — and £n+p+i = 1- The domain (0, 1) of definition of 
the spline functions is the one-dimensional parametric domain. On it, the knot vector H induces 
a partition of the interval (0, 1) that we denote by M. Precisely, we define M as the set of the 
knot spans (£i,Ci+i), i = [p/2] + 1, ■■■,n + [p/2\, that can also be empty due to knot multiplicity 
greater than 1. Empty intervals still play a role in the definition of B-splines and are graphically 
represented as points close one to the other, as proposed in [33J. Note that in this representation 
of M, the number of lines is the knot multiplicity with one exception: for each boundary knot (at 
or 1) of an open knot vector in M we represent only a multiplicity of \_p/2\ + 1, which is (p+l)/2 
lines if p is odd, and p/2 + 1 lines if p is even (see Figure [TJ. The reason for this construction of 
M will be motivated in the next section. 

Finally, it is worth noting the relationship between the space S P (E) and the space of derivatives 
Sp-i(E/), and their respective meshes M and M'. The meshes M and M' may differ only as regards 
the number of points at the boundary. Indeed, according to the definition above, if p is odd both 
meshes coincide, and if p is even the number of elements of M' with respect to M is reduced by 
two, one on each side. 
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3.1.2. Anchors and Greville sites 

In this section we present the concept of anchors and of Greville sites as points in the parametric 
space (0, 1) which may be associated to each B-spline function. Greville sites, which are also known 
as knot averages, are classical and can be found for instance in [30], while the concept of anchors 
has been introduced recently in [52] . 

Since splines are not interpolatory, the association of functions to points (or, as we will see, 
other geometric entities) is somehow more arbitrary than with Lagrangian finite elements. Anchors 
and Greville sites are two different choices, and we present here both. Anchors are very useful 
when dealing with non-tensor product extensions of splines as T-meshes, while Greville sites (and 
related geometric entities) carry 

degrees of freedom in a more natural way. 

Given a B-spline function N i}P ((), and its local knot vector E itP = . . . we set: 

if p is odd, the anchor A associated with Ni tP (() is the central knot of 5j iP . If p is even, the anchor 
A associated with N itP (() is chosen to be the midpoint of the central knot span of namely: 

( A := ^ t+ p/ 2 + &+p/2+i ^ rjn^ p OS j^j on Q f anchors f or degrees p = 2 and p = 3 are represented 



Note that obviously the correspondence between anchors and B-splines functions is one to one, 
but different anchors A ^ A may lie at the same position ( A = ( A . A remedy to this abuse of 
notation, at the cost of more complex definition, is proposed in [UJ where the use of both an index 
and a parametric domain is proposed. 

The set of anchors is denoted as A p = A P (S). When p is odd anchors are located at all knots of 
the partition M (which may be repeated), while when p is even anchors are located at midpoints of 
all elements in M (including the ones of zero area). Indeed, this fact is the reason for the definition 
of M in particular as regards to the multiplicity of boundary knots. 

Most often, we will use anchors to index functions and local knot vectors. Namely, for an 
anchor A G A p , E A and B A (() will denote the corresponding local knot vector and B-spline 
function, respectively. When no confusion occurs, the subindex may be removed. 

Remark 3.1. The B-splines are, in general, not interpolatory at the anchor A 6 A P {E?), while 
they are interpolatory at knots having multiplicity p. This always happens at Q = and ( = I, and 
happens in the interior of the parametric domain where the basis is C° continuous, i.e., at knots 
with multiplicity p. See e.g., Figure^left) . 

Given A e A p , and E A = ..,£ i+p+ i} for some i, then the Greville site is defined as: 



Unlike anchor positions, Greville sites are all different one from the other, when the multiplicities 
irij verify rrij < p and thus B-splines are all continuous. The Greville sites induce a partition of 
the unit interval, referred as Greville mesh and denoted M^- These concepts are ubiquitous in 
spline theory and geometry representation. Greville sites are intimately related to control points 
and control polygon whose properties we briefly recall in the next section. 

3.1.3. B-spline curves 

A B-spline curve V in M. 3 is defined by a parametrization in the interval (0, 1), in the form 



2 



in Figure [T] 




6+i + ... + £ i+p 
p 



(8) 




(9) 
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Figure 2: B-spline curve and its control polygon (left), and the same curve after one step of /i-refinement (right). 



where C A e IR 3 are called the control points. Control points are in a one-to-one correspondence 
with B-spline basis functions. The piecewise linear interpolation of the control points gives the 
control polygon Tc- See Figure [2] for an example. 

The control points C A have an important role not only in the definition of the spline parametriza- 



tion ( 15 ), but also in the visualization and interaction with spline geometries within CAD softwares. 
Indeed, it is common in CAD softwares to represent, together with the parametrized curve T, the 
control points C A and the associated control polygon T C - Typically, the CAD user defines or 
interacts with the control points in order to input and modify the geometry. Since the B-splines 



are not in general interpolatory (recall Remark 3.1), then the control polygon Tc differs from T, 



but it is "close" to it. Precisely, Tc converges to T under /i-refinement. This convergence is proved, 
e.g., in [30] and discussed here below. 

We introduce the usual Lagrangian basis for piecewise linear polynomials on the Greville mesh 
M G , denoted by X A (-): 

The control polygon Tc is then parametrized by the mapping : [0, 1] — > Tc defined by 

Fo(C) = cAaA (0, < £ < 1, (11) 

that is, Fc and F share the same control points. When F is smooth enough, the following 
approximation estimate holds (see, e.g., [301 Ch. XI]): 

sup ||F(C)-F C (C)|| ^h 2 , (12) 
Ce[o,i] 

h denoting the mesh-size. In other words, Tc approximates T up to an error 0(h 2 ) under h- 
refinement. A graphical representation of this convergence can be seen in Figure [2] 

3.2. Multivariate B-splines 

Multivariate B-splines are defined from univariate B-splines by tensor product, see for instance 
P, SO] • Anchors are defined in a similar way. We give here a quick overview. 
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3.2.1. Knot vectors, B- spline functions, anchors, Greville sites 

Let d be the space dimensions (in practical cases, d — 2, 3). Assume ni G N, the degree pe EN 
and the pe-open knot vector Ei = {£e,i, ■ ■ ■ > C^+j^+i} are given, for £ = l,...,cL These knot 
vectors define a tensor product mesh M in the parametric domain Q = (0, l) d where, as in Section 
3.1.2 we have to take into account the knot multiplicity. The multiplicity of a knot vector in 



is represented graphically by lines (d = 2) or planes (d = 3) one close to the other, while the 
boundary is treated exactly as in one dimension. 

The set of anchors is defined on M as the Cartesian product A Plt _ jPd (Ei, . . . , E d ) = A Pl (E 1 ) x 
... x A Pd {Ed). Considering, for example, the trivariate case (d = 3) and recalling the definitions 
from Section |3. 1.2 for the univariate case, we have that: if all pe are odd the anchors lie at the 



vertices of the mesh; if both pi and P2 are odd and p^, is even, then the anchors are middle- 
points of the vertical edges of M; if both pi and P2 are even and p% is odd, then the anchors 
are centers of the horizontal faces of M; if all pi are even the anchors lie at the center of the 
elements of M, and so on. As in the univariate case, the anchors may be located at the center 
of zero length edges or zero area faces or empty elements, according to knot repetition. Also, the 
computation of the local knot vectors for each anchor follows from the univariate case. Given 
an anchor A = (Ai, A 2 , A 3 ) G A pltP2)P3 = A PljP2)P3 (E 1 ,E 2 ,E 3 ), we have that its coordinates are 
(CiS (2, (3) = (C Al ,C A2 ,C A3 )- The three local knot vectors (one in each coordinate direction) 
corresponding to A are defined as Ef := [Ei] Ai for % = 1,2,3 and the B-spline associated to A is 
constructed by tensor product as: 

< P2 , P3 (0 = < 1 (Ci)^ 2 (C 2 )^t(C3)- (13) 
withC = (Ci,C2,C 3 )efi = (0,i) 3 . 



The B-spline functions (13) span the space S Pl!P2!P3 (Ei,E 2 ,E 3 ) (or simply S PliP2>P3 ), which is 
the space of piecewise polynomials of degree pi in the xg direction on M, whose continuity at the 
internal mesh plane Q = £1% is C Pe ~ mi < k , being the multiplicity of in the knot vector Eg. 



To each anchor A e A pltP2jP3 (E 1 ,E 2 ,E 3 ) (or, equivalently, to each B-spline function (13)) we 
also associate a Greville site in the natural way, that is 

7 A =(7i A ,7 2 4 ,73 A ) (14) 

where each is defined as in ([8]), from the local knot vector Ef. Connecting adjacent Greville 
sites, we obtain the Greville mesh Mg, which is a regular tensor product mesh with all elements 
of positive volume. 

3.2.2. Spline and NURBS geometries, multi-patch domains 

Analogously to spline curves, a trivariate single-patch spline parametrization of the domain 
Q C R 3 is F : VI Q defined as a linear combination of B-splines, 

F(C)= c X, P2 , P a(0, wither, (15) 

where C A G M 3 are called control points. In a similar way, it is also possible to define bivariate 
spline domains in IR 2 or surfaces in IR 3 , which are commonly used in CAD (see, e.g., [HI 129]). 

The control points C A have again the same important role in the visualization and interaction 
with geometries within CAD softwares. Now, the concept of control polygon generalizes to the 
control mesh Mc, that is, the mesh connecting the control points. Figure [3] shows an example 
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geometry, with its control points and control mesh. The control mesh defines a polyhedral domain, 
denoted Qc, which is an approximation of Q: again, the control domain Qc converges to Q under 
/i-refinement. This is stated as in the univariate case: we introduce the usual Lagrangian basis for 



Figure 3: Representation of a geometry (green), with its control points (blue) and control mesh (black lines) for 
splines of degree 3. 

piecewise trilinear polynomials on the tridimensional Greville mesh Mg, still denoted by \ A (-), for 
the sake of brevity, 

' 1 if A = A\ 



X A ( 7 A ' 



if A ^ A'. 



The control mesh is the image of the Greville mesh Mg through the piecewise trilinear mapping 
F c : -> Q c , 

f c (c)= c A A A (fl, with c en, (16) 

A^Ap 



which is a parametrization of Qq, since 

,A\ r^A 



Fc(7 A ) = C 



When F is smooth enough, as for (12), we have 

8up||F(0-F c (0ll ^h 2 . (17) 

The control mesh plays a fundamental role in structural mechanics applications where the 
unknowns are sought as displacements of control points. In our work, we will show how this 
interpretation can be used also in other contexts. 

Remark 3.2. When pi = P2 = P3 = 1 (and all anchors have multiplicity one) the Greville sites 
coincide with the anchor representations, i.e., *y A = C, A , and F(£) = F C (C)> VC e Q, that is, Vtc 
and Q coincide. 



9 



In CAD and isogeometric analysis the geometry is often parametrized by Non Uniform Rational 
B-splines (NURBS). NURBS are generated from projective transformations of splines (see [6]). A 
trivariate single-patch NURBS parametrization of the domain Q C M. 3 is a function F : Q — > Q 
defined as quotient of linear combination of B-splines, 

Yja&a C a w a B£ 1P2 (Q 

J?(Q — fc P1-P2>P3 Pl.P2.P3 V»/ £ £ Q (18) 

where C A are the NURBS control points and w A the positive NURBS weights. 

In order to enhance flexibility and allow for more complex geometries, the definition of tensor- 
product spline and NURBS parametrized domain can be generalized to domains that are union of 
N images of cubes; precisely 

closure (fi) = closure (19) 

k=l,...,N 

where the Qk — Ffc(O) axe referred to as patches, and are assumed to be disjoint. Each patch has its 
own parametrization F&, defined on its own spline or NURBS space. The whole Q is then referred 
to as a multi-patch domain. For the construction of discrete fields on a multi-patch domain Q we 



will introduce in Section |4.4| suitable conformity assumptions. These will restrict the framework 
to configurations where it is easy to implement the proper continuity of the fields at the patches 
interface. 

In this paper, Q is assumed to be parametrized either by spline or NURBS functions but the 
unknown fields are always constructed by splines. This means that, in case of NURBS geometries, 
we leave the isoparametric concept which is a fundamental assumption for isogeometric methods 
in the context of continuum mechanics (see [8]). 

4. The spline complex 

This section is devoted to present the spline spaces that are compatible with the De Rham 
complex. The definition of the spaces is taken from [22], and is given in three dimensions (though 
the same construction is generalizable to arbitrary dimension). We first recall the construction on 
the parametric domain Q, and then the discrete spaces on the physical domain Q are obtained by 
the push-forward mapping associated to Q . As shown in [22] , it is also possible to complement 
these spaces with commuting and continuous projectors, in the setting of the so called Finite 
Element Exterior Calculus (see [2]), however this issue is not discussed here. Instead, we discuss 
the selection of a suitable basis for the implementation of the proposed spaces, and the meaning 
of the associated degrees-of-freedom. We will see that the proposed spline spaces are a natural 
high-order extension of classical low-order Nedelec hexahedral finite elements of the first family 
(see [IB]), obtained in this setting for degree p\ = p 2 = P3 = 1, and that in a natural way they are 
related to cochain or chain complexes of the mesh where they are defined. 

4-1. Complex on the parametric domain 



We recall the following property of univariate splines, from Section 3.1.1 the derivative of a 
(continuous) spline is a spline, and in particular 

S P (E) ViP'), (20) 



10 



where E' is the (p — l)-open knot vector that coincides with the p-open knot vector E except for 
the boundary knot repetitions. Moreover, we have that the derivative of the B-spline associated 
to an anchor A in A P (E) is a linear combination of the B-splines associated to the previous and 
next adjacent anchors A~ and A + in A p -i(E') (only one adjacent anchor for the first and last 
A G A P (E))\ precisely, denoting by E^ the local knot vector (formed by p + 2 knots) of A and by 
E p _ x the local knot vectors (formed by p+ 1 knots) of A ± 1 and with the general notation of Section 

|^K-) = |^f^'i](0- 

„_!i are the length of the support of the (p — l)-degree B-splines iVpE^J, that is, the 
difference of the last and first knots in the local knot vectors E^_ x 

holds with iVfS^] 



3.1.1 we have 



V 



(21) 



where |H^_ li 

An example is given in Figure 

When A is the first (resp., last) anchor, ( [21 

a well known property of B-splines (see [6j [30] ) that also suggests the following scaling of the basis 
functions of S P (E) and S p -i(E' 



S P (E) = spaa {!#(.) = N[Ef}(.) : A G A P (E)} 



S. 



p-i\ 



V 



-p-il 



span < D^-) = j^NlE^i.) : A e A 



(22) 
(23) 



The scaling in (23) gives the Curry-Schoenberg B-splines (see [301 Ch. IX]), that have been already 
used in isogeometric analysis in [2 



Indeed, with the bases (22) and (23), the matrix associated 



to the operator ^ is the edge-vertex incidence matrix related to the mesh M, when p is odd, or 
the vertex-edge incidence matrix related to M, when p > 2 is even. We recall that M also contains 
zero length edges and repeated vertices. 




Figure 4: Derivative of the spline associated to the anchor A as a linear combination of the splines associated to 
A~ and A+. 



The observations above are the key ingredients of the trivariate construction. Following [22J, 



and using the notation of Section 3.2, we introduce the following discrete spaces on the parametric 
domain Q 
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xl 



^Pl-l,P2,Piip'11 "2J X "Spi,P2-l,P3 ("lj '-'2 J '-'3) X '5pi,P2,P3-l('-'l> "2; "3)) 

^Pi^-l.Pa-ll"!' "2J "3) X ^Pl-l^Ps-lO^l, "2, —3) x Spi-lj^-lj^G^l, n 2> ^3), 



(24) 



<Spi-i 



Pi-l,P2-l,P3-ll"l5 "2) "3/1 



From (20), they form a De Rham complex: 



0. 



(25) 



Moreover, we have the following result. 



Theorem 4.1. The sequence (25) is exact. 

Proof. This result has been already presented in [22]- We present an alternative proof, that will 
be useful Section 15.31 



We have to show that in (25) it holds 

K = ker( grad ), 

im( grad ) = ker( curl 
im( curl ) = ker( div ), 
im(div) = X%. 



(26) 

(27) 
(28) 
(29) 



In particular, we have to prove the inclusion D in (26)-(29), since the other inclusion C is trivial 
in all cases. It is also trivial that 



D ker( grad 



Let u = (ui,U2,U3) G Xl, then define 



Ci 



<2 



C3 



0(0,(2,(3)=/ £1(77,0,0)^77+/ u 2 ((i,rj,0)dv+ u 3 (( u (2,v) d V 
Jo Jo Jo 

it is easy to check that u = grad</> when curlu = 0, and that (j) e then 

im( grad ) D ker( curl ). 
Consider tp e X%, and define v = (tJi, 0, 0) G such that 

ui(Ci) C2, Ca) = / ¥{11,(2,(3) drj, 
Jo 

as before, it is easy to check that = divv and that v e X\\ then 

im(drv) D 



(30) 



(31) 
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In order to complete the proof we need to show that 

im( curl ) D ker( div ), 

which is implied by 

dim(im( curl )) = dim(ker( div )). 



To count dimensions recall from Section 3.1.1 that dim(S Pe (Ei)) = tig, dim(S Pe -i(E' e )) 
then from Section 3.2 and from (24) we get 



(32) 
ni - 1; 



dim(X°) 
dim(^) 
dim(X, 2 ) 
dim(X 3 ) 

Then, by ph-p^. 



= nin 2 n 3 , 

= (m - l)n 2 n 3 + m(n 2 - l)n 3 + nin 2 (n 3 - 1), 

= ni(n 2 - l)(n 3 - 1) + (m - l)n 2 (n 3 - 1) + (ni - l)(n 2 - l)n 3 , 

= ( ni -l)(n 2 -l)(n 3 -l). 

dim(im( curl )) = dim(X^) — dim(ker( curl )) 
= dim(X^) — dim(im( grad )) 
= dim(X^) - dim(X°) + dim(R) 
= 2n 1 n 2 n 3 - n 2 n 3 - n i n 3 - n>in 2 + 1 



(33) 



and by (29) 



dim(ker( div )) = dim(X^) — dim(im( div )) 
= d\m(X 2 h ) - dim(X, 3 ) 

= 2{m - l)(n 2 - l)(n 3 - 1) + (n 2 - l)(n 3 - 1) + (m - l)(n 3 - 1) + (m - l)(n 2 - 1) 
= 2nin 2 n 3 - n 2 n 3 - riin 3 - n x n 2 + 1, 



which gives (32), and as a consequence (28). 

We now show how suitable basis functions for the spaces can be constructed and as well 
associated to geometric entities of the mesh M by using the concept of anchors. We focus on basis 
functions first, and inspired by (22)-(23) we define them as follows: 



X° h = span{C i y b£(Ci)B£(C2)B£(C S ) : A = (A U A 2 ,A 3 ) G A Pl , P2 , P3 (E l7 E 2 ,E 3 )} , (34) 

Xl = span IUIIU III, with 

I = {C H- ^-iCCO^CCaJ^CCsJei : A = (A 1 ,A 2 ,A 3 ) G ^-^(Si, H 2 , S 3 )} , 
// = {C H- ^(CO^-iCCO^CCs)^ : ^ = (i4i.i42.i4a) G ^^(Sx, 5 2 , 5 3 )} , 
77/ = {C H- ^(Ci)^ , (C2)^ , -i(Cs)3s : ^ = (Ax.Aa,^) G ^^^(Si, ~ 2 , S' 3 )} , 

= span 7 U 77 U 777, with 

7 = {C ^ B^CO^-ifo^-ifo)* : A = (A 1 ,A 2 ,A 3 ) G ^Lpi.pa-i.ps-i^i, H 2 , H 3 )} , 
77 = {C ^ ^_i(Ci)^ , (Ca)^-i(Cs)3 2 : A = {A X ,A 2 ,A 3 ) G y^-i^-i^, E 2 , E 3 )} , 
777 = {C ^ < 1 _ 1 (Ci)^ 2 _ 1 (C 2 )^ 3 3 (C3)e 3 : A = (A 1 ,A 2 ,A 3 ) G y^-i^-i^i, E' 2 , E 3 )} , 
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Xl = span{C H- I^_i(C0^-i(C 2 )i^-i(C8) : A = (A U A 2 ,A 3 ) e ^-i^-i^-^Hl, S^, H^)} , 

(37) 

where {e£}^ = i j2) 3 denote the canonical basis of M 3 . We remark that all basis functions defined in 



(34)- (37) are non-negative. 



We discuss now the association of the anchors of the bases (34)-(37) to the mesh M that is 



associated to X°, that is, obtained from the knot vectors Hi, H2, H 3 . We focus on the relevant case 
p = pi = p 2 = P3 and consider two possible choices: p is odd, or p is even. 

When p is odd, as an immediate consequence of the definition of anchors in one space dimension, 
we have that: 

• anchors associated with X® are A PiPtP (Ei, H 2 , H 3 ), which are located at the vertices of M, i.e., 
there is one degree of freedom per vertex; 

• anchors associated with X\ are located at edges of M and there is one degree of freedom per 
edge. Indeed, e.g., anchors associated with the first component of the space X\, which is 
5p-i,p,p(S'i, H 2 , H 3 ), are A p -\ >p #(z! x , H 2 , H 3 ) and are located at the edges oriented as This 
means that to each edge of the mesh a is associated a basis function tangential to the edge. 

• anchors associated with X\ are located at faces and there is one anchor per face. More in 
detail, if we consider the first component of X%, which is S , p iP _i iP _i(Hi, H' 2 , H' 3 ), the corre- 
sponding anchors are A PiP -i )P -x(Ei, H 2 , H 3 ) and are located at the barycenter of the faces / 
such that / is orthogonal to ei. This means that a basis functions normal to the face is 
associated to the face. 

• anchors associated with X\ are yLp_i iP _i >p _i(H / 1 , H' 2 , H 3 ) and located at barycentres of all 
elements of M. 

We now turn to the case when of even de gree y > 2, p\ = p 2 = p 3 = p. Note that, according 



to our definition, and as explained in Section 3.1.1, the meshes corresponding to the spaces X^, 
X\ and X\ differ from M due to the different number of repeated lines at the boundary. Instead 
of working with different meshes for different spaces, equivalently, we represent in this case too 
the anchors of all spaces on the mesh M of keeping into account only the interior geometrical 
entities for the representations of anchors of X\, X\ and X\. 

Using the definition of anchors we immediately deduce the following: 

• anchors associated with X® are at the barycentres of all elements in M; 

• anchors associated with X\ are attached to the barycentres of internal faces of M and the 
corresponding vector basis function is normal to the face; 

• anchors associated with X\ are attached to internal edges of M and the corresponding vector 
basis function is tangent to its corresponding edge; 

• anchors associated with X\ are attached to internal vertices of M. 

Clearly, the positivity of the bases induces an orientation of edges and faces of the mesh M. 



With the bases (34)-(37), the discrete differential operators in (25) are represented by incidence 
matrices for the corresponding geometrical entities. If p = p\ = p 2 = p 3 is odd, then the operator 
grad is represented by the edge-vertex incidence matrix of M and when p > 2 is even, by the face- 
element incidence matrix of M. We observe that, unlike in compatible finite elements, the matrices 
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representing the differential operators in the selected bases (34)-(37) are essentially independent 
of the degree. 

The fundamental consequence of the observations above is stated in the following proposition. 
Proposition 4.2. The following holds: 



The spline complex (25) for odd degree p is isomorphic to the cochain complex associated 
with the partition M. 



The spline complex (25) for even degree p is isomorphic to the chain complex associated with 
the partition M without its boundary, that is, when only the interior geometrical entities 
(faces, edges and vertices) are taken into account, as seen above. 



As a matter of fact, this observation, together with the structure of the matrix representation 
of differential operators, makes the geometry of the spline complex for odd degree p very similar, 
if not equal, to the one of the finite element complex of lowest order. However the spline complex 
for p > 1 delivers an approximation which is far superior than the one of low order finite element. 

For p even we have instead a chain complex without explicitly constructing the dual mesh, 
which has no analogue in the finite element framework. 

Moreover, the use of anchors and the structure of the mesh at the boundary guarantee that 
in both the chain and cochain complex the boundary is treated in a simple and canonical way. 
In the case of finite elements this is not case (see e.g. [SOI EI], or (53] ) and, moreover, 

these features can hardly be obtained in conjunction with high-order finite element techniques. 
Discretization methods based on the use of both chain and cochain complexes in the framework 
of isogeometric methods are very promising and object of on-going research. 

We conclude this section with a remark on boundary conditions. Consider the case when 
homogeneous boundary conditions are imposed on the whole boundary dQ, leading to the definition 
of the discrete spaces Xj h := X Q h n iZ£(fi) , X^ h := X l h n H (curl; fi) , X 2 j?t := X 2 n H (div; fi) 

and X$ h :=X%. These spaces are constructed as usual, by removing the functions with non-null 
trace at the boundary, because univariate B-spline functions are interpolatory at the boundary, as 



we have discussed in Remark 3.1 The associated De Rham complex 



is exact, as easily follows by a variation of the argument of Theorem AA_ The same holds in more 
general cases, for example when the boundary conditions are imposed on a part of dfl formed by 
the union of some faces of the cube Q. Since boundary conditions do not represent a conceptual 
difficulty, in order to keep the presentation as clear as possible often in our presentation we will 
not take them into the framework. 

4-2. Push- forward to the single-patch physical domain fl 



Following Section 3.2.2, we suppose that the domain Q is obtained from Q through a spline or 



NURBS single-patch mapping F. Clearly, we need to choose the space for F. 

Assumption 4.3 (Isogeometric mapping). We assume that F is either a spline function in 



[X®] 3 , or F is a NURBS function as in (18), with numerator in [X°] 3 and weight denominator in 
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Assumption 4.3 is indeed very natural in the context of isogeometric methods: it means that the 
discrete fields are constructed from the geometry knot vectors and bases, possibly after refinement. 

We denote by M the image of M through the mapping F. M is then a partition of the physical 
domain Q, similar to the finite element mesh, even though it contains elements of zero area due to 
knot multiplicity. 



h , . . . , Xl on the physical domain Q can be defined from the spaces ( 24 ) 



The discrete spaces X° 

on the parametric domain Q by push-forward, that is, the inverse of the transformations defined 
in Q, that commute with the differential operators (as given by the diagrams @ and (|3 



+ x° h 



grad 



curl 



curl 



> xl 



div 



div 



> xl 



> xl 



that is, the discrete spaces in the physical domain are defined as 



-+ 



(39) 



-> 0, 



xl 


■■= {0 


PW) e x° h }, 


xl 


= {u 




xl 


:={v 


^ 2 (v) e xl}, 


xl 


= w 


G X, 3 }. 



(40) 



We remark that the space X\, which is a discretization of H(curl; Q) , is defined through the curl 
conserving transformation t 1 , and that the space X\, which is a discretization of H(div; Q) , is 
defined through the divergence conforming transformation l 2 . These are equivalent to the curl and 
divergence preserving transformations that are used to define edge and face elements, respectively 
(see [HI Sect. 3.9]). 

Thanks to the properties of the operators ([T]) the push-forwarded spaces X®, . . . , Xl inherit the 
same fundamental properties of X®, . . . , Xl, that we have discussed in the previous section: 



they form an exact De Rham complex without boundary conditions 



x2 ^ xl 



or with boundary conditions 



curl 



> xl 



div 



> xl 



-> o, 



(41) 



xS h 



grad j curl v2 

y A o,/i y A 0,/i 



div 



^ Xqk 



(42) 



the basis functions for X®, 



, Xl are defined by push-forward of the basis functions of 
Xl, . . . , Xl, similarly to (40), and are in one-to one relation with the images of the anchors 
through F. See Figure [5j 



since (39), the matrices associated with the differential operators grad, curl and div on 

f2 are the same as the matrices of grad , curl and div on tt, that is, incidence matrices of 
the mesh M. 



when p is odd (even, respectively), the complex (X®, . . . ,Xl) is isomorphic to the cochain 
(chain, respectively) complex associated to the partition M. 



16 



Figure 5: We show the mesh M and the image of the anchors related to the space X° on an example geometry. 



Finally, the discrete spaces X®, . . . , X^ inherit from their pull-back X®, . . . , X^ optimal approx- 



imation properties, if the geometrical mapping F satisfies Assumption 4J3 and its inverse is smooth 
enough (see [22] for details). 

4-3. Control fields and degrees- of-freedom interpretation 

In this section, we introduce the concept of control fields, thanks to which we give an interpre- 



tation of the degrees-of-freedom of the isogeometric fields defined in Section |4.1| - [4.2| The control 
fields are for the B-spline isogeometric fields what the control mesh is for the the B-spline geometry. 
We recall that from the geometry control points we define Fc (see (16)), the piecewise trilinear 
function on the Greville mesh Mg. The image of Fc is the so-called control domain flc- The 
so-called control mesh Mc (which is a partition of Qc) is the image through Fc of the Greville 
mesh Mg. 



As described in Section 3.2.2 the standard way to manipulate a spline parametrization F is by 
moving its control points, that is, the vertices of the control mesh Mc- The parametrization Fc or, 
equivalently, the control mesh Mc, carries the degrees-of-freedom for the geometry. The distance 



between the two parametrizations F and Fc is at most 0(h 2 ), as in (17). We now apply the same 
rationale for the complex (X®, X%). Let us first focus on scalar functions on the parametric 
domain Q, i.e., on the space X®. Given a spline 



?«)= E c X p ,p(0, withCe^, 



(43) 



where c A are its control variables, we associate the piecewise trilinear function defined on the mesh 

5o(C)= E cAxA (0, wither, (44) 



which carries the same degrees-of-freedom for cf) and indeed is close to (the distance between 
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the two functions is at most 0(h 2 ), analogously to (12)). By this relation, we can interpret the 
degrees-of-freedom c A of as the values of 0c at each Greville site in M^. 

It should be noted that, if the values of these degrees-of-freedom are chosen wisely, splines 
deliver approximation error of order 0(h p ) in the norm of H 1 (Q), where p is the degree of the 
spline, while the corresponding trilinear function can only provide approximation errors of order 
0(h). 



Let now the geometry come into play Using (40), we set: 



oF 



and 



5C 



oF 



c 



?C- 



(45) 



The degrees-of-freedom for are the values of 0c at the vertices of Mc, that is, at the control 
points. Or, we can say that the field 0c determines, or controls, 0, and its degrees-of-freedom 
are the values of 0c at control points. In Figure [3j the location of control points (blue bullet) 
is shown on an example geometry. The field 0c plays the role of control field for 0. As for the 
parametric space, there are wise choices of the degrees-of-freedom which ensure an approximation 
error of order 0(h p ) in the norm of H 1 (Q), while the corresponding trilinear function can only 
provide approximation errors of order 0(h). 

The same reasoning can be applied to the whole complex (X®, ...,X^), defined in Section |4.1 
and |4.2[ from degrees vt. and knot vectors H^. Indeed, we introduce the control complex (Z t 

and 4.2, with the following choices for Z®: 



from degrees pi and knot vectors_ 
which is obtained, still following Section 



4.1 



70 



degrees in all directions equal to 1; 

the knot vector in the £ direction is the ordered collection of points {•yf : A £ A Pe (EA}, 
i = 1,2, 3, along with repeated and 1 to make the knot vectors open, 



and replacing the geometric mapping F with Fc in the pullbacks Q. The complex (Z®, ■■■,Zf k ) 
corresponds to the low order finite element complex defined on the control mesh Mc and it is 
immediate to see that if m (45) belongs to X®, then the corresponding <pc belongs to We 
denote by 7° : X° — > Z\ the operator which associates to 0c, and in an analogous way, we define 
the operators P h : X 3 h — > Z J h , j = 0, ... ,3. These operators are represented by identity matrices 
when the spaces are endowed with the bases described in Section |4~T] It is not difficult to see that, 



in view of the structure of the matrices associated to differential operators, the following diagram 
commutes: 



-T 



grad 



curl 



4 



xl 



div 



-T 



xl 



-> 



4 



(46) 



7° 



grad 



curl 



7 2 



div 



7 3 



0. 



Let us comment about the meaning of the diagram (46). First of all, it says that the geometric 
structure of the spline complex (X®, ...,X^) is the one of the low order finite element complex 
(Z®, Zfy on the control mesh. The discrete fields in X ] h can be associated to control fields in 
Z 3 h , through the operator P h , as we have discussed for j = above. For example, we can say 
that there is a Nedelec field % which controls u and the degrees-of-freedom are, in this case, its 
circulation on the ed ges of the control mesh Mc. Moreover, following a reasoning similar to the 
one in Section 3.2.2 the operators P h are point-wise converging to the identity when h goes to 
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Figure 6: Representation of the degrees-of- freedom location for the space X^, on the green geometry for degrees 
Pi = Vi = P3 = 3. 

zero. We stress again that the order of approximation of the complex (X°, ...,X^) is 0{h p ) while 
the control complex Z|) only exhibits first order convergence in the norm of X\ 

Finally, it should be noted that, as it is well known, the complex .., Z%) is always isomorphic 



to the cochain complex of the partition Mc, while for the complex (X®, ..,X%) Proposition 4.2 
holds. This is in accordance with the fact that, when pi are all even, the control mesh Mq can be 
interpreted as a partition dual to M, in the sense that the chain of M is isomorphic to the cochain 
complex of M^. 

4-4- Push-forward to the multi-patch physical domain Q 

In this section we construct the spline complex on a multi-patch geometry by addressing the 
questions of how conformity is imposed at the interfaces between patches. 

We consider now a multi-patch domain Q which is parametrized from a reference patch Q 



through the spline or NURBS mappings k = 1,...,N, as in (19). Each patch is endowed 



with a (possibly different) spline space and therefore for each k = 1, . . . , N we can define discrete 



spaces [X°]fc, . . . , [X^]fc such that a De Rham complex (25) holds. Assuming each verifies 



Assumption |4.3[ then, as shown in Section 4.2 we push-forward patch-by-patch the discrete spaces 



[A°] fc , . . . , [Xf] k and obtain, on each Q k = F&(fi), the discrete spaces [A°] fc , . . . , [Xf\ k that fulfill 



the De Rham complex (4.2) on each patch. 



Then, the last and main step is to assemble the spaces X 3 h C @ fe=1 ^\X-h\ki and add the 
relevant continuity conditions at the inter-patches boundaries: trace continuity for X°, tangential 
trace continuity for X^, normal trace continuity for X%, no continuity for X%. For this purpose, 
we introduce a conformity condition as, e.g., in This condition guarantees that the geometry 



parametrizations of the patches are equivalent at the patch interfaces and, since Assumption 4^ 
it can be stated on the spaces [A°]fc. 

Assumption 4.4 (Geometrical conformity). On each non-empty patch interface F = dQ k fl 
dVL k t , with k ^ k! , the spaces [A°] fc | r and [X°] fc /| r coincide, as the corresponding bases do. 
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This means that the meshes and match on T, and therefore 



M= [J M f 



k=l,...,N 



is a locally structured but globally unstructured mesh M on Q. In a similar way, the patch control 
meshes [Mcj^ match conformally and 



M c = |J [M ( 



'CJfc 

fc=l,...,JV 



is a locally structured but globally unstructured mesh Mc of hexahedra on f2c, the union of the 
patch control domains [Ocjfc. 



Assumption 4.4 corresponds to the full-matching conditions of [26], to which we refer for further 
details. 

Having conformity we can implement the continuity conditions easily. Indeed, due to the 



definitions in Sections 4.1-4.3, the needed continuity holds for (X°, ■ ■■,X^) if and only if it holds 
for (Z%, Zfr) on the global mesh Mc. Continuity for Z%, Z\ and Z\ is imposed by merging the 
coincident degrees-of-freedom at the interfaces, which in the case of Z\ and Z\ also requires to 
take into account the orientation; see Figure [7j This is however the same as in finite elements 



(indeed, the control fields are classical finite elements). This merging automatically gives the 



degrees-of-freedom for fields in (X®, ...,Xf 



5. Beyond the tensor product structure: T-splines 

In this section, we generalize the definition of tensor-product B-splines to T-splines [7J [311 [32] . 
The theory of T-splines is well developed in two dimensions (see the very recent papers [371 GB1 
W7\ [54] ) while it is still incomplete in three dimensions (though some recent important advances 



have been recently proposed in [IH]). For this reason, we only present in Sections 5.1 and 5.2 



T-splines in two dimensions and construct, in Section 5.3, a discrete T-spline based complex. The 



extension to three dimensions is given in Section 5.4 by tensor-product of the two-dimensional 



T-spline spaces with B-spline one-dimensional spaces. 
5.1. T-mesh 

Let ne G N and the degree pi G N, and let = {&,i, . . . , £t,ni+p t +i} be a p^-open knot vector for 
£=1,2. A T-mesh is a rectangular tiling of the unit square [0, l] 2 , such that all vertices belong to 
^! x S 2 . A T-mesh may contain interior vertices that connect only three edges, called T-junctions, 



that break the tensor product structure of the mesh (see Figure 10). We will say that a T-junction 
is horizontal (respectively, vertical) if the missing edge is horizontal (resp. vertical). By an abuse 
of notation, we still denote a T-mesh by M. 

As before, we represent the knot multiplicities by repeated lines close to each other, with now 
the line multiplicity possibly varying along lines (see [321 Section 4.3]). The only exception are the 
boundary lines, that maintain the same multiplicity all along the line. As in B-spline meshes, the 
vertical (resp. horizontal) lines at the boundaries have multiplicity |_pi/2j + 1 (resp. |j>2/2j + 1). 
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(a) Control variables of the two patches before merging. (b) Control variables of the merged patch. 

Figure 7: Implementing continuity for X\ on a two-patch domain. The orientation of the edges at the interface is 
chosen as that of the lower patch. 

5.2. Analysis suitable T-meshes and T-splines. 

We define, for a horizontal (resp. vertical) T-junction T, the k-b&y face-extension as the 
horizontal (resp. vertical) closed segment that extends from T in the direction of the missing 
edge until it intersects k lines of the mesh M. The A;-bay edge-extension is defined analogously 
extending the segment in the opposite direction. 

Following [17], we define the extension of a horizontal (resp. vertical) T-junction T the union 
of the LG°i + l)/2j-bay face-extension, and the [(pi — l)/2j-bay edge-extension (resp. the union 
of the [(p 2 + l)/2j-bay face-extension, and the [(p 2 — l)/2j-bay edge-extension). More precisely, 
if pe is odd we extend (pe + l)/2 bays in the direction of the missing edge, and (pp — l)/2 bays in 
the opposite direction; if pp is even we extend pp/2 bays in both directions. An example is given 
in Figure [8j 

Definition 5.1. A T-mesh M is analysis suitable for degrees p\ and p 2 , denoted M G AS Pl)P2 , if 
vertical extensions and horizontal extensions do not intersect. 

Analysis suitable T-meshes were first identified in [37] in the bi-cubic case, and generalized 
to arbitrary degree in [UJ. Despite their very geometric definition, analysis suitable T-meshes 
and T-splines enjoy fundamental properties which make their use in isogeometric analysis really 
promising. Some of these properties will be discussed in what follows. 
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Figure 8: Extensions for degree p\ 



= 2 (horizontal) and P2 = 3 (vertical). Dashed lines represent the face extensions. 



1/6 2/6 3/6 4/6 5/6 1 



Figure 9: Representation of the local knot vectors for degrees p\ = 2 and P2 — 3. The local knot vectors are 
3f = {0, 0, 1/6, 2/6}, Ef 1 = {0, 0, 0, 1/6, 2/6}, and Sf = {3/6, 4/6, 5/6, 1}, Ef = {0, 2/6, 3/6, 4/6, 5/6}. 



As in the case of B-splines, anchors are inferred from the T-mesh M and their position depends 
upon the parity of p\ and p^- For example, if both pi and pi are odd, anchors are at the vertices 
of the M, if they are even, anchors are at the barycenters of elements and so on (see [E]). We will 
denote the set of anchors by A PltP2 (M), or simply A PltP2 . 

T-spline basis functions are constructed as B-splines associated to the anchors A Pl ^ P2 (M.), and 
defined from two local knot vectors, Ef = {£ 1)ix , . . . ,£i,i P1+2 } C Hi and Ef = {£ 2)il , . . . ,£ 2 ,i P2+2 } C 
H2. To construct the horizontal knot vector Ef we trace the horizontal line through A and select 
its intersections with vertical lines of M, depending on the degree p\. if p\ is even we choose the 
first (pi + 2)/2 intersections to the left of A, and the first (pi + 2)/2 to the right; if p\ is odd we 
first select the coordinate of the anchor A, and then the first (p\ + l)/2 intersections to the left 
and to the right of A. In the case that we arrive at the boundary, we add the value or 1 as many 
times as needed to complete the p\ + 2 entries of Ef. The construction of Ef is analogous, and 
depends on P2. Examples are shown in Figure [9] for p\ = 2 and P2 = 3. For more details we refer 
to [321. 
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The T-spline function associated to the anchor A is denoted as: 



Bi 



Pl,P2 



(o = NiEfmNimm, c = (Ci,c 2 ) e (o,iy 



(47) 



they are linearly independent (see [37]) and by definition span the T-spline space T pltP2 = T Pl P2 (M): 



T piiP2 (M) := span{£? 



A 

P1,P2 



AeA puP2 (M)}. 



(48 



Definition 5.1 guarantees fundamental properties of the T-spline space (48). In [38j I37J it is 



defined a dual basis for the T-spline functions constructed from an analysis suitable T-mesh, thus 



proving the linear independence of (47) (see also [37]) and good approximation properties for the 
space (48). We also remark that the construction of the local knot vectors described above is 



analogous to the one in jUJ for analysis suitable T-meshes. 

Finally, we define the extended T-mesh of M, and denote it by M ext , as the T-mesh obtained 
from M by adding all the T-junction extensions. The extended T-mesh, sometimes also called 
Bezier mesh, is the minimal mesh such that the functions (47) restricted to the non-empty el- 



ements are bivariate polynomials of degree (pi,P2)- The importance of the extended mesh for 
implementation |55j, local refinement [53] and approximation properties of T-splines |47j is al- 
ready known. In particular, for the implementation, and in order to ensure accuracy, integration 
has to be performed on the elements of M ext and this means that the data structure is constructed 
based on M ext . 

Finally, a key result which is useful in the construction of compatible T-spline discretizations, 
is the characterization stated in the following proposition. 

Proposition 5.2. Given an analysis suitable T-mesh M, if furthermore no T-junction extensions 
of any kind intersect each other or intersect mesh lines with multiplicity greater than one, then the 
T-spline space (48) is the space of all piece-wise bivariate polynomials of degree (pi,P2) on M ex t 



with the same continuity of the T-spline functions (47) at the mesh lines 



Proof. The case p\ = pi = 3 has been covered in [54J, while the general case is a work in progress 
by A. Bressan in [SB]. Related works are also [ST], and [SB] which show the mathematical complexity 
of the problem. The condition that the extensions do not intersect lines with multiplicity greater 
than one can be removed at the price of a more complex statement, which we do not consider here 
for the sake of simplicity. 



5.3. Two-dimensional De Rham complex with T-splines on the parametric domain (0, l) 2 

The aim of this section is to introduce a two-dimensional T-spline based De Rham complex, 
thus generalizing the tensor-product construction of section |4j Throughout this Section we will 
assume, for the sake of simplicity, p\ = P2 = p. The results are also valid in the general case, but 
the proofs become more intricate. 

As for B-splines, T-splines spaces are constructed by a suitable selection of the polinomial 
degree in the two directions and by a suitable design of the mesh, that is, the knot vectors. The 
main difference now is that we need to modify the mesh M, depending of the form degree, not 
only at the boundary but also around T-junctions. 

Let p G N, let Si,H2 be two p-open knot vectors, and let M G AS PiP be a T-mesh with knot 
repetitions, as defined in Section 5.1 The starting mesh is M° = M, on which we define the space 
of scalar fields: 

Y°:=T p , p (M°). (49) 
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The T-splines vector fields are defined on the two T-meshes 7v[\ and Mjj. If p is odd, M.\ 
is obtained from M by adding the first-bay face-extension of all horizontal T-junctions. If p is 
even, is equal to M everywhere but on the boundary where, due to the definition of M, and 



recalling Section |3.1.1[ the first and the last vertical columns of elements of M are removed. We 
define analogously M-j,, reasoning in the vertical direction: if p is odd M2 is defined by adding the 
first-bay face-extension of all the vertical T-junctions, and if p is even, it is defined by removing 
the first and last horizontal rows of elements of M. Then the vector fields and the rotated vector 
fields are defined as 

:= T p _ 1)P (M|) x T PtP ^{M\). (50) 

?£* :=T p ^ x {M\)xT p _ 1}P {M\). (51) 

Finally, the last space is defined on the T-meshes M 2 : if p is odd, M 2 is obtained from M 
by adding all the first-bay face-extensions (horizontal and vertical), and if p is even it is defined 
by removing the first and last rows and columns of elements in M. 

Then 

Y h '■= T p _i iP _i(M 2 ). (52) 



An example of the sequence of meshes is shown in Figure 10 for p = 3, and in Figure 11 for 



p = 2. We notice that whenever M° = M is a tensor product mesh, the construction is equivalent 



to the one presented in Sect ion |4T] for B-splines. Indeed, for odd p the four meshes are equal to M, 
because there are no T-junctions, and for even p they only differ in the number of line repetitions 
on the boundary. 

The choice of these meshes becomes clear when computing the derivatives. For instance, let 
p = 3 and consider the simple example of a mesh with only one horizontal T-junction, as in 



Figure 12(a) Choosing the anchor A G A PtP (M°) located at the T-junction, it is clear fro m ([21~| 



that - is a linear combination of B^_ x and B^_ x , with A\ A r G yL p _ ljP (Mj) as in Figure 12(b) 
Hence, ^ G T p _ hp {M\) (see d50b), but since A 1 ^A p . hp {M°), we have £ T p _i iP (M°). The 



argument is analogous for the partial derivative with respect to the y direction, with vertical 
T-junctions. 

We have the following result. 

Proposition 5.3. Assuming M° G AS P)P , it holds that M} G AS p _i )P , 3Vkj> G AS PjP _i ; and M 2 G 
ASp-i.p-i. Moreover, (M°) cxt = (M\) cxt = (Ml) ext = (M 2 ) cxt . 



Proof. The result is an immediate consequence of M G AS PiP , and the length of the extensions 



specified in Section 5.2 



Remark 5.4. Note that, although the four meshes are different, all integral computations are 
carried out in the extended T-mesh, which is the same for all the spaces. As a consequence the 
four spaces can be implemented within the same data structure, which is based on one single mesh, 
but with different basis functions for each space. This is also what occurs with standard finite 
elements. 
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(a) Mesh M° 
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(b) Mesh 



(c) Mesh Ml 



(d) Mesh M 2 

^ure 10: Sequence of meshes for the spline complex, with their respective anchors, for p 
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(a) Mesh M° 



(b) Mesh 



(c) Mesh Ml 



(d) Mesh M 2 

Figure 11: Sequence of meshes for the spline complex, with their respective anchors, for p = 2. 
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(a) Anchor A E A Pt p(M°) 



(b) Anchors A l ,A r £ A p -i, p (M\) 



OB 



Figure 12: The first partial derivative d F x F is a linear combination of B^_ t p and B p _ x p . 



The bases for Y®, . . . Y£, are formed by T-spline functions (47) with a scaling as in Section 4.1 
Precisely, introducing the notation D[Ef](Q) := jJ^N[Ef](Q), we have: 



K, = Bpan{(Ci,Ca) iV[~f](G)iV[~^(C 2 ) : A G A P , P (M )} 



(53) 



Yfr = span / U //, with 

/ = {(&,00 ^ D[S^](Ci)iV[H^](C2)ei : A G ^^(M})} , 
J/ = {(Ci,C 2 ) ^ iV[H^](Ci)P[H^](C 2 )e 2 : A G A p ^{^\)} . 



Y,l* = span / U II, with 

/ = {(Ci,Ca) ^ ^[Hi 1 ](Ci)^[^](C2)e 1 : A G . 
// = {(Ci,Ca) ^ ^[Hf](Ci)iV[^](C 2 )e 2 : A G -Vi^Mj)} . 

F, 2 = span{(Ci,C 2 ) ^ £>[^](Ci)£>[H^](C 2 ) : A G *A p -i, P -i(M 2 )} 
The main result of this section is the following. 



(54) 

(55) 
(56) 



Theorem 5.5. Under the assumptions of Proposition 5.2 , the following two-dimensional com 
plexes 

rot 



r> grad ^ 



(57) 

_ h r* h * k , u, (58) 

where rotu = (<9iw 2 — <9 2 Wi) i/ie scalar rotor and rotw = (d 2 u, —diu) T is the vector rotor, are 
well defined and exact. 



-> 0, 



-* 0, 



Proof. In the proof we only consider (57), since (58) is equivalent. The well posedness of the 
complex follows from 

and ^oT : 9£ ->■ Y^, (59) 



;rad : Y® -> F, 1 



which, in turn, easily follows from the definitions (49)-(52) and from Proposition 5.2 
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Exactness of (57) means 



K = ker( grad ), 
im( grad ) = ker( rot ), 
im( rot ) = Y%. 



(60) 
(61) 
(62) 



The first part, i.e., (60), is obvious. Moreover (61) is also simple: indeed if u e has null rot 



then u = grad0, where, e.g., 

0(Ci,C 2 ) 



ux(r),0)dr) + / u 2 ((i,r)) drj. 

JO 



(63) 



Since, u = grad <p, then has to be element by element (of M, ext ) a p-degree tensor-product 
polynomial. Then, <p inhe rits the interelement regularity from u and has the one of functions in 
Y®. Then, by Proposition 5.2, <p e Y®. The last point, (62), follows from the dimension formula 

dim(i; ) + dim(Y£) = dim(F h 1 ) + 1. (64) 



Indeed, using (61), (60), and (64) 



dim(im( rot )) = dim(Y h 1 ) — dim(ker( rot )) 

= dim(y /i L ) — dim(im( grad ))) 

= dim(Y^) — dim(yj°) + dim(ker( grad ))) 

= dim(i > h L ) - dim(y ft °) + 1 

= dim(? /t 2 ). 



In order to prove (64), we recall the Euler's formula for the T-mesh M 

F + V = E + 1, 



(65) 



where Fq is the number of faces, Eq the number of edges and Vo the number of vertices of M, 
including knot repetitions, zero length edges and empty elements. The proof is different for odd 
and even p. 

Case 1) Let p be odd. We can separate the edges into horizontal and vertical ones, and with self- 
explaining notation we have E = Eq + Eq . Similarly, the vertices can be divided into horizontal 
T-junctions, vertical T-junctions and all the other vertices (including those on the boundary), in 
the form Vo = V H + V V + V Q + . For odd p the meshes M.\ and M2 are constructed by adding 
the first-bay face-extension of horizontal and vertical T-junctions, respectively. Thus, using the 
assumption that T-junction extensions do not intersect, the number of horizontal edges in Mj is 
E^ = Eq + V Q , and the number of vertical edges in M2 is E\ = Eq + V V . Similarly, the mesh 
M 2 is contructed by adding all the first-bay face-extensions, and the number of faces in M 2 is equal 
to F 2 = F + Vq H + V V . 



Since p is odd, and from the positions of the anchors in every mesh (see Figure 10), the 
dimensions of the spaces are 



dim(F /l °) = Vo, dim(F, 1 ) = E? + E\ 



v 



E q + Vq H + Vq V , dimft 2 ) 



^2 = ^0 + Vq H + V c 



> 



and using (65) the proof is finished. 
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Case 2) Let p be even. We denote by V B and E B the number of boundary vertices and boundary 
edges in M, and we note that V B = Eq . As before, we distinguish between horizontal and vertical 

V FT R B B H 

edges, Eq = Eq + Eq , and also for the boundary edges Eq — E ' + E ' . For even p the mesh 
(resp. M2) is constructed by removing the first and last columns (resp. rows) of elements from 
M. Hence, the number of vertical edges in M} is E\ = Eq — Eq ,v , and the number of horizontal 
edges in M2 is E B = Eq — E ' . Similarly, the mesh M 2 is constructed by removing the first and 
last rows and columns of elements from M, thus the number of vertices in M 2 is V2 = Vq — V B . 
From the position of the anchors for even p (see Figure [TT|) , the dimensions of the spaces are 

B 
• 



dim(y°) = F , dim(r, 1 ) = E\ + E? = E - E B , dim^ 2 ) = V 2 = V - V, 



Using (65) and that V = Eq the proof is finished. □ 



5-4- Three-dimensional De Rham complex based on T-splines and B-splines 

We construct a three-dimensional complex on the parametric domain by tensor product of the 



two-dimensional T-spline complexes (57)-(58) and the one- dimensional complex (20). Then we 
define the spaces on the parametric domain Q = (0, l) 3 : 



xl 
xl 
xl 



: [Y£ ® S P (E)] x [Y° ® S P ^(E% 



(66) 



which form a complex of the kind ^ (or ^ if we also impose homogeneous Dirichlet boundary 
conditions). 

Assume now th at t he geometry map F is tensor-product single-patch spline or NURBS, and 
fulfills Assumption 4.3[ now with X° defined as in (66). Therefore, the push-forwards ([2])-([3]) give 
the correct complex (X®, . . . , X%) on Q: this procedure is completely analogous to what we have 



already described in Section ^2 and is not detailed here. 

It is not a difficulty to consider, more generally, a multi-patch, or a T-spline geometry mapping. 
This is not detailed here, for the sake of brevity, but the first case will be addressed in the numerical 
tests of the next section. 

5.5. Concluding remarks on the T-spline complex 

As it appears from our presentation, the understanding of the T-spline complex is much less 
sound than the one of the spline complex, even in two space dimensions. Moreover some of the 
properties we have studied for splines do not hold in general for T-splines. For example, 

• the matrices corresponding to the operators are no more the incidence matrices of the mesh 
M and a similar fact is true for standard finite elements with hanging nodes, i.e., the T-spline 
complex with p = 1; 

• the definition of control mesh and control fields is not trivial especially when p is even and 



the analogue of Section 4.3 is not available for T-splines. This deserves further studies. 
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6. Numerical results 



In this section we present numerical tests showing the behavior of isogeometric methods for 
electromagnetic problems. Since numerical tests for B-splines have already been presented in other 
works, see e.g., [59| 127 ] 122]. we will concentrate here on examples involving also T-splines. All our 
numerical tests have been performed with the Matlab library GeoPDEs [60]. It should be said 
though that GeoPDEs does not have full T-splines capability, and in particular does not provide 
any T-splines adaptivity in the sense of [53] . 



The mappings we use in this section always verify the Assumption in Section 5.4 , and can 
be either single-patch or multi-patch; the meshes we describe are the ones corresponding to the 
space X%. The meshes for the other spaces are constructed following the procedure detailed in 



Section 5.3 In the figures of this section, repeated lines of the mesh are represented with thicker 
lines, independently of the number of repetitions. In all cases, internal mesh lines have multiplicity 
one. 

6.1. Maxwell eigenproblem in the square domain 

As a first test we solve the two-dimensional eigenvalue problem: Find (u, u) G H (rot; Q) x R 
such that 

/ rot u rot v = u 2 / u ■ v Vv G H (rot; Q) , (67) 
Jq Jq 

in the square domain Q = (0, it) 2 , for which the exact eigenvalues are u 2 = m 2 + n 2 , with m,n = 
0, 1, . . .. The aim of this test is to show that the discretization of the problem with T-splines does 
not present spurious modes. 

The coarsest mesh consists of 8 square (non-empty) elements in the left half, and 4 rectangular 
(non-empty) elements in the right half, thus creating several T-junctions on the vertical line £i = 



0.5. Finer meshes are created by dividing each element into 4 (see Figure 13). 

In Table [T] we present the first non-null eigenvalues for degree 3 and for the sequence of meshes 
explained above. The results show that there are no spurious eigenvalues, and that a good con- 



vergence rate is obtained. In Figure 14 we display the first non-null eigenvalues computed with 
discretizations of degree 4 and 5 in a mesh formed by 768 non-empty elements, and its comparison 
with the exact eigenvalues. Again, it is seen that the discrete eigenvalues are computed with the 
right multiplicity. 



Figure 13: Coarsest mesh for the square, and mesh after one refinement step. 
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□ Degree 4 
O Degree 5 
Exact 
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Rank in the (non-null) spectrum 



Figure 14: First non-null eigenvalues computed in the square for degrees 4 and 5. 



Mode 


Exact 


Computed 


(1,0) 


1.00000 


1.00001 


1.00000 


1.00000 


1.00000 


1.00000 


(0,1) 


1.00000 


1.00005 


1.00000 


1.00000 


1.00000 


1.00000 


(1,1) 


2.00000 


2.00016 


2.00000 


2.00000 


2.00000 


2.00000 


(2,0) 


4.00000 


4.00396 


4.00004 


4.00000 


4.00000 


4.00000 


(0,2) 


4.00000 


4.03882 


4.00134 


4.00002 


4.00000 


4.00000 


(2,1) 


5.00000 


5.00395 


5.00003 


5.00000 


5.00000 


5.00000 


(1,2) 


5.00000 


5.10164 


5.00208 


5.00002 


5.00000 


5.00000 


(2,2) 


8.00000 


8.05454 


7.99989 


8.00001 


8.00000 


8.00000 


(3,0) 


9.00000 


9.06255 


9.00135 


9.00001 


9.00000 


9.00000 


(0,3) 


9.00000 


9.12399 


9.02102 


9.00057 


9.00001 


9.00000 


(3,1) 


10.0000 


10.0614 


10.0014 


10.0000 


10.0000 


10.0000 


(1,3) 


10.0000 


10.2361 


10.0324 


10.0007 


10.0000 


10.0000 


(3,2) 


13.0000 


12.8159 


13.0028 


13.0000 


13.0000 


13.0000 


(2,3) 


13.0000 


13.2002 


13.0091 


13.0004 


13.0000 


13.0000 


(4,0) 


16.0000 


17.9413 


16.0181 


16.0002 


16.0000 


16.0000 


(0,4) 


16.0000 


19.8934 


16.2962 


16.0076 


16.0001 


16.0000 


(4,1) 


17.0000 


19.9586 


17.0181 


17.0002 


17.0000 


17.0000 


(1,4) 


17.0000 


20.8937 


18.0245 


17.0092 


17.0001 


17.0000 


(3,3) 


18.0000 


21.4707 


18.7373 


18.0008 


18.0000 


18.0000 


(4,2) 


20.0000 


24.0689 


20.0191 


20.0002 


20.0000 


20.0000 


(2,4) 


20.0000 


26.1844 


21.6138 


20.0056 


20.0001 


20.0000 


d 


oi. 


74 


184 


548 


1852 


6764 


number of zeros 


21 


65 


225 


833 


3201 



Table 1: First non-null eigenvalues computed in the square for p = 3. 




Remark 6.1. We have also solved the previous test by the mixed formulations in [61], that make 
use of the full two-dimensional De Rham complex (57). The computed non-null eigenvalues are 
the same as for the plain formulation (67), while the zero eigenvalues are filtered with the mixed 
formulation. These results, that we do not present here for the sake of brevity, confirm that the 
construction of the De Rham complex with T-splines is correct. 
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6.2. Maxwell eigenproblem in the thick L-shaped domain 

As a second test case, we solve the three-dimensional eigenvalue problem: Find (u, oS) G 
H (curl; fl) xR such that 

I curlu- curlv = cj 2 / u v Vv G H (curl;fi), (68) 
Jq Jq 

in the thick L-shaped domain Q — E x (0, 1), where £ = (—1, l) 2 \ [—1, 0] 2 . From |B2], it is known 
that the reentrant edge introduces a singularity in the first eigenfunction, which only belongs to 
the space H 2 / 3 ~ e (Q) for any e > 0. 

It is well known that in order to recover the optimal convergence rate we need to suitably refine 
the mesh toward the reentrant edge, see e.g., [63J or [64]. Anisotropic elements need to be used 
in this case (see [21] for some theoretical background on the topic). We propose here a dyadic 
refinement based on T-splines. 

For the geometry representation, the thick L-shaped domain is parametrized as the union of 



three cubic patches. Following Section 4.4, scalar fields in X® are only continuous at the interfaces 



between patches, and the fields in X^, which are used in the discretization of (68), are only 
tangentially continuous at these interfaces (like for standard edge finite elements), but at least 
C p ~ 2 within patches. 

The refinement is obtained via T-splines by dyadic partitioning of elements which are close to 
the reentrant edges [651 Ch. 4]. We perform the refinement first in an L-shaped two-dimensional 
section, and then propagate to the three-dimensional domain with a uniform mesh in the z- 



direction, as already explained in Section 5.4 The refinement is performed identically for every 
patch, in such a way that conformity can be ensured at the patch interfaces. 

To construct the two-dimensional mesh, at each refinement step, and for every patch, we refine 
a small square region near the reentrant edge, subdividing each element into 4. Then some T- 
junction extensions are added, depending on the degree, to make the mesh analysis suitable, as 



defined in Section |5.2| For instance, in the example of Figure [15] we start with an 8 x 8 mesh for 
each patch, which is drawn in black. At the first step we refine a region of 3 x 3 elements on each 
patch. Since the degree is p = 4, two-bay extensions must be added to make the mesh analysis 
suitable. The refined elements at this step are given by the blue lines. At the second refinement 
step, which is marked in red, we first refine a square region of 2 x 2 elements, and again we add the 
two-bay extensions to make the mesh analysis suitable. Finally we remark that, since the dyadic 
partition and the definition of analysis suitable T-mesh depend on the degree, different meshes are 
used for different degrees. 

The problem has been solved for degrees 4 and 5. In Tables [2] and [3] we present the first 
non-null computed eigenvalues in the three cases, and its comparison with the exact solution. In 



Figure [16] we display the convergence rate for the first eigenvalue, comparing the results obtained 
with T-splines and with a B-spline discretization of the same degree in the corresponding refined 
tensor-product dyadic mesh. As can be seen, with T-splines we obtain the same error with an 
important reduction in the number of the degrees of freedom. 
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Figure 15: Example of a two-dimensional mesh in the L-shaped domain £, and its extension to the three-dimensional 
domain Q = £ x (0, 1) , for p = 4. 




10 s 10 4 10 5 10 s 10* 10 s 



Degrees of freedom Degrees of freedom 

(a) Degree 4 (b) Degree 5 



Figure 16: Convergence of the first eigenvalue in the thick L-shaped domain. 



Exact 


Computed 


9.63972384472 


9.64482260735 


9.64055367165 


9.63986647533 


9.63977706731 


9.63974511214 


9.63972731966 


11.3452262252 


11.3444193267 


11.3450973393 


11.3452056015 


11.3452178503 


11.3452228875 


11.3452256921 


13.4036357679 


13.4036330719 


13.4036359208 


13.4036359870 


13.4036357431 


13.4036357654 


13.4036357699 


15.1972519265 


15.1973643408 


15.1973310163 


15.1973301300 


15.1972556440 


15.1972524223 


15.1972523662 


19.5093282458 


19.5144198480 


19.5101576732 


19.5094708082 


19.5093814308 


19.5093494993 


19.5093317180 


19.7392088022 


19.7392474090 


19.7392464705 


19.7392464473 


19.7392098765 


19.7392090606 


19.7392090522 


19.7392088022 


19.7392474115 


19.7392464714 


19.7392464480 


19.7392098765 


19.7392090617 


19.7392090536 


19.7392088022 


19.7392854949 


19.7392835833 


19.7392835402 


19.7392109156 


19.7392092780 


19.7392092574 


21.2590837990 


21.2591164396 


21.2591199815 


21.2591200740 


21.2590848357 


21.2590840605 


21.2590840611 


d.o.f. 


4042 


7126 


11018 


16162 


22630 


34894 



Table 2: First non-null eigenvalues computed in the thick L-shaped domain for p — 4. 
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19.7392095886 


19.7392095763 


19.7392095758 


19.7392095758 


19.7392088273 


19.7392088022 


19.7392095886 


19.7392095763 


19.7392095758 


19.7392095758 


19.7392088273 


19.7392088022 


19.7392103678 


19.7392103456 


19.7392103444 


19.7392103444 


19.7392088514 


21.2590837990 


21.2590824582 


21.2590845213 


21.2590845754 


21.2590845767 


21.2590838179 


d.o.f. 


5891 


9883 


14827 


20723 


28105 



Table 3: First non-null eigenvalues computed in the thick L-shaped domain for p = 5. 



6.3. Maxwell source problem in three quarters of the cylinder 

For the third test case we consider the model source problem: Find u G (curl; fl) such 

that 

/ curlu- curlv+ / u • v = / f ■ v Vv G H ^ D (curl; fl) , (69) 
Jn Jn Jn 

where Ho,r D (curl; fl) is the set of functions with null tangential trace on C <9fl, i.e., 
Ho,r D (curl; fl) := {v G H(curl; fl) : v x n = on r D }. 



The geometry fl is three quarters of a cylinder of radius and length equal to one (see Figure 17), 
that in cylindrical coordinates is given by fl = {(r, 9,z) : < r < 1, < 9 < |7r, < z < 1}. 
We impose the null tangential component on = {(r, 9,z) : 9 G {0, |vr}}, and the source term 
f is taken such that the exact solution is u = grad (r 2 / 3 sin(2#/3) sin(7r^)), i.e., it is singular in 
the first two directions, but it is regular in the z direction, for which the local refinement of the 
previous test is well suited for this case. 

As in the previous example, the domain is defined with three patches, and the discrete fields 
X\ C H(curl; fl) are only tangentially continuous between them. The construction of the mesh in 
the parametric domain is identical to the one in the previous example: for each patch we first create 
a two-dimensional mesh locally refined towards the corner, and extend it to the three dimensional 
domain by tensor product. The mesh is then mapped to the physical domain, as can be seen in 



Figure 17 



The problem is solved with T-splines of degree 3, and also with B-splines of the same degree 
in the corresponding tensor product mesh. The errors in H(curl)-norm for the two methods are 



compared in Figure 18 As in the previous example, with T-splines we are able to obtain results 



similar to those of B-splines with a reduction of the number of degrees of freedom. 

6-4- Numerical simulation of a twisted waveguide 

As the last numerical test we use T-splines to simulate the propagation of a singular mode in 



a waveguide with a twist. The configuration, which is presented in Figure 19(a) , is the same given 
in [66l Ch. 8], changing the material discontinuity by a geometric inhomogeneity (the twist). The 
problem is solved in a waveguide with a twist of 90 degrees, with a section of three quarters of 
the circle of 2 cm radius, and for which we assume that the walls are perfect electrical conductors. 
We also assume that the waveguide extends to infinity without other inhomogeneities, and it is 
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Figure 17: Example of a mesh for three quarters of the cylinder. 




Degrees of freedom 

Figure 18: Comparison of the error for T-splines and B-splines. 



truncated by the planes Ti and T 2 to obtain the computational domain, which consists of three 



different regions: a middle region where the waveguide is twisted (see Figure 19(a)), and two 
straight regions near the ports, to keep the inhomogeneity far enough from them, in such a way 
that only the dominant mode TEio can propagate without attenuation. The total length of the 
computational domain is 24 cm: 4 cm for each straight region, and 16 cm for the region with the 
twist. The frequency uj is taken equal to 32 GHz, and it is between the cutoff frequencies for the 
first mode (21 GHz) and the second mode (33.84 GHz). 

Following [66] , and working in the time harmonic regime at a given frequency u, the (complex- 
valued) electric field E G H ,r D (curl; Q) is solution of the problem 

[ ( curl E ■ curl G - k 2 E ■ G) + f ift (n x E) • (n x G) = 

2i/3 10 / E inc -G, VGe Ho,r D (curl;fi), 

where k = ^u 2 HqEq with /i and e the magnetic permeability and electric permittivity of free 
space. The incident electric field E mc at the port F±, and the wavenumber of the first mode /3io 
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are defined as 



E iDC (x,y,z) = ei (x,y)e 



10 



1.2 _ 1.2 



In the case of waveguides of rectangular or circular section, the value of the constant kio and the 
mode eio are known. In the general case, they can be obtained by solving a 2D eigenvalue problem 
on the port T\, which consists on finding the minimum eigenvalue fcio € K, and its associated 
eigenvector e 10 G H ( rot ; Ti), such that 



rot e 10 rot v 



ft io 



ri 



eio-v Vv G H (rot;ri). 



(71) 



The electric field E in equation (70) is discretized with T-splines of degree 3, using the approach 



already explained in Section 5.4 The two-dimensional T-mesh for the section is built as in the 
previous examples, and in the z direction we use one element for each straight region near the 
ports, and 4 elements along the twist, for a total of 7936 degrees of freedom. For the solution of 



the 2D problem (71), it is enough to restrict a field in H(curl;fl) to its tangential components 



on the port r\, which in practice is equivalent to solve with two-dimensional T-splines. 



The magnitude of the real part of the computed solution E is shown in Figure 19(b), which 



shows that the mode is correctly propagated. Finally, we also compute the reflection and trans- 
mission coefficients, given by the equations 



-i/3z 



R 



eio 



-2i/3zi 



/r, 



T 



eio 



-2i/3z 2 



eio ■ eio 



/r, 



eio ■ eio 



and we obtain the values \R\ = 0.0025 and |T| = 
affect the propagation of the mode, as expected. 



0.9998, which confirms that the twist does not 




22,332 

(a) Geometry of the waveguide (b) Real part of the computed electric field 

Figure 19: Geometry of the waveguide and real part of the computed electric field. 
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